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Abstract. - We present a lattice-based numerical method to describe the non equilibrium be- 
havior of a simple fluid under non-uniform spatial conditions. The evolution equation for the 
one-particle phase-space distribution function is derived starting from a microscopic description 
of the system. It involves a series of approximations which are similar to those employed in theo- 
ries of inhomogeneous fluids, such as Density Functional theory. Among the merits of the present 
approach: the possibility to determine the equation of state of the model, the transport coef- 
ficients and to provide an efficient method of numerical solution under non-uniform conditions. 
The algorithm is tested in a particular non equilibrium situation, namely the steady flow of a 
hard-sphere fluid across a narrow slit. Pronounced non-hydrodynamic oscillations in the density 
and velocity profiles are found. 



Introduction. — The behavior of a confined fluid can 
be different from that of a bulk fluid in many important 
Aspects. First of all the confinement induces density inho- 
rnogeneities which may determine a variety of phenomena 
leaving no counterpart in bulk systems [1]. The presence 
of surfaces, not only alters the average equilibrium prop- 
erties of fiuids, but also affects their time-dependent be- 
havior such as diffusion, momentum and energy transport, 
j^s a result, a great effort is currently devoted to the un- 
derstanding of fluid physics at the nanoscale (see [2] and 
rleferences therein). 

In the last thirty years, a massive effort has been de- 
voted to the understanding of system at thermodynamic 
equilibrium and new techniques have been developed, 
among these Density Functional theory (DFT) being per- 
haps the most versatile [3]. On the contrary, we do not 
have a similar control over the behavior of non equilib- 
rium systems. Dynamical extensions of DFT have been 
introduced and tested successfully in the case of colloidal 
suspensions, where the damped character of the dynamics 
makes the density the only relevant field [4]. Instead, in 
the case of standard fiuids one needs to fully account for 
the momentum and energy transport. The natural pro- 
cedure seems therefore to consider the evolution of both 
positions and momenta of the molecules and to derive an 



equation for the phase-space distribution able to convey 
the structural information about the microscopic nature of 
the fluid. Such a task can be achieved by using a modified 
Enskog-Boltzmann approach, which has been proposed by 
different authors [5-9]. 

Our present goal is to provide a procedure able to de- 
scribe simultaneously the discrete nature of fluids and the 
transport properties at interfaces and in confining geome- 
tries. To this purpose, we briefly recall that the Lattice 
Boltzmann (LB) method represents a very eflacient scheme 
derived from kinetic theory to simulate fluid flows and 
transport phenomena [10,12]. Being based on the contin- 
uum Boltzmann equation, it accounts for a faithful repre- 
sentation of the macroscopic hydrodynamic behavior. The 
original idea of LB is to model fiuid flows by simplified 
kinetic equations, the Bhatnagar-Gross-Krook (BGK) re- 
laxation operator being a popular choice to simulate the 
Navier-Stokes evolution at macroscopic level [10,11]. The 
LB of McNamara and Zanetti [13] averages the micro- 
dynamics by solving the kinetic equation of the particle 
distribution instead of tracking the motion of each par- 
ticle. Whereas the original formulation was designed to 
describe lattice gas of particles and the attention was fo- 
cused on large scale properties, modern versions of the 
LB aim to incorporate a more realistic behavior of the 
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fluid, fn particular, since the pioneering work of Shan and 
Chen [f 4, f 5] there has been a large effort to include the 
effect of microscopic interactions between fluid molecules 
into the description [6, 7, 9, 16, f 7]. Such interactions have 
been accounted for by considering the intermolecular force 
field at a particular point in the fluid. This modification is 
sufficient to describe non-ideal gas behavior such as phase 
coexistence and phase separation and various surface phe- 
nomena, but it is still unsatisfactory since it docs not allow 
for a systematic prediction of the macroscopic properties 
starting from the microscopic level. 

In other approaches, in order to circumvent such a dif- 
ficulty the collision terms have been dealt with explicitly, 
but at the price of introducting gradients of the macro- 
scopic fields [6,7,9]. 

Inspired by recent work on inhomogeneous fluids [18,19] , 
we propose an alternative scheme to evaluate the interac- 
tion terms involved in the kinetic equation. As a result, 
we obtain evolution equations for the distribution function 
whose structure is very similar to that employed in Den- 
sity Functional theory, without the explicit evaluation of 
gradients of macroscopic fields. We believe that the strat- 
egy of using coarse-grained quantities instead of gradients 
of the relevant fields may give a better representation of 
fluids confined to narrow systems, as we have learned from 
the study of equilibrium fluids, where gradient expansions 
usually give poorer results than non-local density func- 
tionals. 

In addition, the theory allows to compute equilibrium 
quantities, such as the surface tension, together with ac- 
curate estimates of the bulk transport properties. 

Theory. — In general, the intermolecular potential of 
a simple fluid can be decomposed into a repulsive, respon- 
sible for the microscopic structure at high packing frac- 
tion, and an attractive contribution, which plays a major 
role in determining the thermodynamic properties. How- 
ever, in the present paper, which serves as to introduce the 
method, we confine ourselves to discuss a harshly repulsive 
fluid, the hard sphere system. 

The evolution of the phase-space one-particle distribu- 
tion /(r, V, i) is represented by the following transport 
equation: 

dtfir, V, t) + V V/(r, V, t) + ■ V^„/(r, v, t) - 

m 

m{r,^,t) (I) 

where Fg is an external force and fJ[/] represents the ef- 
fect of the interactions among the fluid particles which we 
describe via the revised Enskog collision operator [20] 

n[f] = a''-^ J dv2 y"dae(a-vi2)(a-vi2) x 

{52(r, r - a, t)/(r, v[,t)f{r -a,v'^,t)- 

52(r, r + a, t)/(r, vi, i)/(r + a, V2, t)} (2) 

where v[ and are scattered velocities determined from 
Vj = vi — (ct • vi2)o' and V2 = V2 + (ff • vi2)o', cr is the 



hard-sphere diameter, a is the unit vector directed from 
one particle to another, and g2(ri,r2|n) is the positional 
pair correlation function. 

Equation (f), together with (2), represents a nonlinear 
evolution equation for the one-particle distribution func- 
tion, /(r,v,t), and could in principle be solved numeri- 
cally by brute force, or by considering a suitable trunca- 
tion of the open hierarchy of equations for the moments of 
the distribution function. However, since one is chiefly 
interested in the evolution of the hydrodynamic mo- 
ments of the distribution function n(r, t) = J dvf{r, v, t), 
n(r, <)u(r,t) = /c?vv/(r,v,t) and ^n{r,t)kBT{r,t) = 
^ / dv(v — u)^/(r,v,t), it is possible to further simplify 
the form of Q without altering the local transfer of momen- 
tum and energy. To this purpose, following Dufty, Santos 
and Brey (DSB) [5], we replace (f) and (2) by a simpler 
equation, where the complicated interaction between non 
hydrodynamic modes is approximated via a BGK-like re- 
laxation term — j/o[/(r, v, i) — n(r, t)(/)A/(v)], vq being a 
phcnomcnological collision frequency, chosen as to repro- 
duce the kinetic contribution to the viscosity. The DSB 
equation can be cast in the form: 

dtf{r,v,t) + vV/(r,v,<) + 5^.V„/(r,v,i) = 

TO 

- t^o[/(r,v,i) - ?i(r,<)(/)M(r,v,i)] 
+ m/?(/)M(r,v,i){(v-u).C(i)(r,t) 

+ (3) 

where /3 = l/kBT{r,t) and 0M(r,v, = 
f . p/^cxpf Moreover 

Cj^^ (r, t) = [ d^rv.n ^--Y. y.P^r, t) (4) 
J m ^ 

where P^^ is the coUisional transfer part of the pressure 
tensor, and 




^ - - E [^'^'^^ ^) + E P^>^ (r , t) V. (r , t)] 

(5) 

arises from the coUisional contribution to the heat flux. 

Equation (3) reproduces the correct form of the hydro- 
dynamic equations, but treats in an approximate fashion, 
viz. within a relaxation time approximation, the higher 
velocity moments contributing to /(r, v, t) and associated 
with kinetic modes [5] . 

The specific forms of C^^-* and C^^^ needed to ren- 
der feasible the method will be given below. Before do- 
ing that, we notice that formula (4) establishes a con- 
nection between the present method and Density Func- 
tional theory. In fact, in the case of vanishing velocity 
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and temperature gradients, eq. (4) can be expressed with 
the help of the excess Helmholtz free energy, Texcin], as 
C(i)(r,<) = --^n(r,tW^pf4r- On the other hand, the 
term C^^''{r,t) describes non isothermal processes and has 
no counterpart within the standard DFT formalism. Un- 
der equilibrium conditions, i.e. in the absence of velocity 
gradients and thermal gradients one can derive the follow- 
ing equation of state with the bulk pressure given by 



Pbulk 



1 



bulk 



rib - 



(6) 

where we set the kinetic contribution to the pressure ten- 
sor P^™ = "1 / dv{v — u)^{v — u)fif{r, V, t) to be diagonal 
and equal to the ideal gas value. In the general non equi- 
librium case, we consider eq. (3) by taking into account 
the full contribution to the pressure tensor stemming, e.g., 
from shearing perturbations. 

Methods and Results. — We shall solve eq. (3) by 
means of the lattice Boltzmann method [12] which has the 
advantage of being very efficient and robust and accomo- 
dates in a natural fashion informations about the micro- 
scopic model. Wc restrict our attention to the isothermal 
situation. 

The original LB method can be viewed as a discretiza- 
tion procedure in velocity space of a kinetic equation 
[6,21,22]. We will employ the popular D3Q19 lattice, 
constituted by a set of 19 populations. The continuous 
velocity v is replaced by a set of discrete velocities c^, 
with 1 = 0, 18, which are vectors which connect neighbor- 
ing mesh points r on a lattice and where the null vector 
Co accounts for particles at rest. Accordingly, the distri- 
bution function is replaced by an array of 19 populations, 
/(r,v,i)^/,(r,t). 

The systematic procedure to discretize the ki- 
netic equation is based on expanding the dis- 
tribution in a finite set of Hermite polynomial 

/(r 



v,0 = u{v)J2!Lo^^\^^t)h'i!M/vT'^l^- where 
(27r4)-2''e-"'/24, vj. = ^kBT/m and by 



considering a Gauss-Hermite quadrature of order 2G to 
evaluate the integrals of the type 



dvv°''...v°"'' /(r,v,t) 



E 

i=0 
G 

i=0 



7(r, C,;,t)/w(Q) 



(7) 



where Wi are quadrature weights and /i(r, i) = 
'Wif{r,Ci^t)/uj{ci). Therefore, the kinetic moments are 
evaluated via 



K 



(8) 



Analogously, we expand the collision operator on the finite 
Hermite set, 

n(r, V, t) = u;{v) f2 i^^f (r, t)hl^ W 



and evaluate its moments via 

C«(r,i)= / dvf7(r,v,i)/i«(v) 



(9) 



The propagation of the populations is achieved via a 
time discretization to first order and a forward Eulcr up- 
date, 

/,(r + c„ t + 1) - /,(r, t)=w,Y, (r, (10) 

Concerning the BGK term appearing in eq. (3), stan- 
dard calculations lead to the following expression for the 
local equilibrium 



n(r, t)$A/(c») 



1 + ^ V'c,j,,Uj,(r, t) 



(11) 



which is tantamount to a low-Mach [0\Ma?'\) expansion. 

In order to obtain a practical scheme we have to spec- 
ify the explicit form of the coUisional contributions to 
the pressure and to the heat flux. This is done by con- 
sidering the fact that in a hard-sphere fluid, momentum 
and energy fluxes can be transferred instantaneously even 
when the velocity distribution function has a Maxwellian 
form, provided it peaks at the local hydrodynamic velocity 
with a spread determined by the local temperature. The 
explicit expression for the coUisional contributions then 
yields cf^ = 0, fulfilhng mass conservation, and 

C^^^(r,t) = -v\ J d^aa„g2{r,r + a\n)n{r,t)n{r + a,t) 

fcB[r(r + a,0-r(r,t)] ] 
+ 2m4 / ^^^^ 

The next C^^) term, governing energy transport, can be 
derived explicitly and reads 

C(^^(r,t) = v^cr'^ J d^ag2{r,r + a\n)n{r,t)n{r+a,t) 



1=0 



1 kB[T{r + a,t)-T{r,t)] ^ 

^TT mvT 1 



(13) 
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Having, now, an explicit representation for Co and 
C'^' in terms of the local values of the five hydrodynamic 
fields we can solve by iteration the equation of evolution 
for /(r, V, i), compute the new values of the fields and 
proceed. Notice that in order to extract the momentum 
flux and the heat flux one has to use formulae (4) and (5) 

The inclusion of energy transport in LB is, however, 
conditioned by the stability of the numerical method. In 
fact, it is well known that non- isothermal LB schemes suf- 
fer from rather severe instabilities already at the level of 
ideal fluids [12]. In the following, we have thus considered 
the regime in which thermal gradients are small, i.e. by 
taking T = const in eqs. (12) and (13). 

The radial distribution function appearing in (12) is 
constructed according to the following prescription, dat- 
ing back to Fischer and Methfcssel [23]. At first, one de- 
fines a coarse-grained density n{r, t) via a uniform smear- 
ing over a sphere of radius cr/2, and the coarse-grained 
packing fraction is fj{r,t) = TTa^n{r,t)/6. Next, the equi- 
librium radial distribution function, g2(j,r + a\n), is re- 
placed by the following approximation [24] (72(1", r + c) ~ 
[(1 - fjir + a 12) 12 + fy2(r + a/2)/4]/[l - + a I2)f . 

In order to evaluate the surface integrals we choose u to 
be an even multiple of the lattice spacing and employ a 18- 
point quadrature over a spherical surface [25]. With this 
choice, the elements arising from gi and the hydrodynamic 
moments are taken from 6 on-lattice quadrature points 
while the elements arising from the remaining 12 off-lattice 
points are constructed via a linear interpolation from the 
surrounding on-lattice elements. 

We have found that this approximation for the pair cor- 
relation function gives excellent results when the coupling 
between density and current in eqs. (12) and (13) is ab- 
sent, corresponding to a dynamical DFT treatment. How- 
ever, when the coupling is active, correlations are spuri- 
ously enhanced even for the static case. This problem can 
be traced back to the well-known compressibility error in- 
trinsic to the LB discretization [12], i.e. the fluid current 
is not rigorously divergence free in absence of external 
forcing. In order to alleviate such a problem, we have 
introduced in the density-current convolution of eq.(12) a 
space-dependent regularizing factor (1 — exp(— |r— r^jl/cr)), 
where is the position of the wall. An alternative to this 
intervention could be the use of an Hermite basis set with 
higher components than the one associated to the D3Q19 
lattice scheme. 

Numerical validation. As a first test of our scheme we 
have determined the shear viscosity of a uniform system 
by performing a linear analysis of the equations of motions 
for the hydrodynamic fields (to be published elsewhere), 
obtained from eq. (3) and eq. (12). It can be shown that 
the present theory gives a shear viscosity identical to that 
predicted by a method proposed long ago by Longuet- 
Higgins and Pople (LHP) [26]. In fig.l we display the 
numerical results together with the predictions of the LHP 
theory and the Enskog theory which, as it is well known. 



gives a value of 77 larger than that obtained from the LHP 
route [27]. 
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Fig. 1: Collisional contribution to the excess part of the 
bulk shear viscosity obtained through: the Longuet-Higgins 
and Pople theory (77 = ^,/k^T^g 2{a){nl ^ikfa*) [26] 

(dashed line), Enskog theory (r, = ^^1/^^(0.8^71^!^ + 

0.77375r2(cr)^^(7ij„jj,)^) [27] (dot-dashed line), and numerical 
data obtained from the time decay of a transverse velocity per- 
turbation (squares). 

Next, we have considered the equilibrium structure of 
the hard-sphere fluid. No-slip boundary conditions are 
enforced on the fluid velocity at the wall via a bounce- 
back method [12]. The corresponding density profiles are 
reported in fig. 2 at various values of the average density 
together with a comparison with the results of equilib- 
rium Monte Carlo simulations. Clearly, the fluid is more 
inhomogeneous at higher densities and displays oscillatory 
behavior in the regions adjacent to the walls. These oscil- 
lations become less evident toward the center of the slit. 
As compared to the exact Monte Carlo solution, the LB 
data provide slightly less correlated profiles. It is worth 
mentioning that the LB solution is achieved with a CPU 
effort about 30 times smaller than needed to generate con- 
verged Monte Carlo data. 

By considering the fluid flow induced by the presence 
of a uniform field, Fe, parallel to the walls of the slit, 
the streaming velocity profiles corresponding to different 
values of the bulk density are shown in fig. 3. We notice 
that as the bulk density increases also the current profiles 
display a non-monotonic structure close to the walls. In 
addition, the average streaming velocity decreases as the 
density increases as a consequence of mutual steric hin- 
drance among particles. 

Finally, in fig. 4 we have considered the dependence 
of the streaming profiles on the width of the slit and, 
for the sake of comparison, displayed the results against 
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Fig. 2: Static density profiles in a slit of width Lja — 5 at 
nj„,j. = 0.256 (lower curve), nj^j^ = 0.512 (mid curve) and 
^buik = 0.609 (upper curve) as compared to Monte Carlo sim- 
ulations. The data at higher density have been shifted by an 
arbitrary value. 

the parabolic velocity profiles d la Poiseuille, predicted by 
the Navier-Stokes theory. Interestingly, the central region 
does not exhibit significant deviations from the quadratic 
behavior of the classical theory, while a non monotonic 
profile next to the walls becomes more pronounced for the 
narrower systems [28] . 

Conclusions. — In summary, wc have presented a the- 
oretical analysis and computational scheme which bridge 
hydrodynamics with microscopic structural theories of fiu- 
ids. We stress that the present method has been derived 
starting from a microscopic model of the fluid and, un- 
like the very popular and succesfuU Shan-Chen method, 
truly represents a bottom-up approach to the description 
of fluid transport. 

In order to simulate microscopic flows, however, the LB 
method should be employed with care. In nanofluidic con- 
ditions, the ratio between the mean free path and the 
characteristic length (the Knudsen number) can be rather 
high. On the contrary, the LB method in the current im- 
plementation (with the D3Q19 finite Hermite basis) is de- 
signed to deal with low-Knudsen conditions, namely in the 
range of validity of the Navier-Stokes equations. Recently, 
however, a generalization of the numerical method to deal 
with high-Knudsen conditions has been proposed [29] , that 
basically extends the Hermite basis up to the third order 
and removes off-basis contributions. We believe that such 
an extension, together with the kinetic approach described 
in this paper, would yield a relevant step forward in the 
study of out-of-equilibrium microscopic flows. 

A straighforward generalization of the present work can 
accomodate attractive forces, which are essential ingredi- 
ents in order to describe multi-phase behavior. A detailed 
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Fig. 3: Poiseuille velocities under the influence of the external 
field Fe/a^ = 2.410"*^ for a slit of width L/a = 10 at various 
bulk densities; nl^if. — 0.064 (circles), 0.128 (squares), 0.256 
(diamonds), 0.384 (triangles up), and 0.512 (triangles left). 

description of the calculations reported in the present com- 
munication will be published elsewhere. 
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